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Modeling  and  Analysis  of  Stochastic  Dynamics  and 
Emergent  Phenomena  in  Swarm  Robotic  Systems  Using 

the  Fokker-Planck  Formalism 


Abstract 

This  research  lays  down  a  mathematical  framework  within  which  a  system  of  interacting  entities,  such  as  a 
swarm  robotic  system,  operating  in  a  dynamic  and  uncertain  environment  can  be  analyzed  for  its  properties  such  as 
stability,  robustness,  and  emergent  behavior.  The  framework,  based  upon  stochastic  differential  equations  and  the 
Fokker-Planck  formalism,  allows  the  calculation  of  transient  and  steady  state  probability  densities  of  the  states  of 
these  systems.  The  framework  has  been  applied  to  two  scenarios  related  to  swarm  systems.  The  first  scenario  relates 
to  a  robotic  swarm,  the  dynamics  of  which  is  inspired  by  ant  foraging.  The  article  develops  a  continuous  time 
reaction-diffusion  model  based  on  Keller-Segel  model  of  bactetrial  chemotaxis  and  establishes  a  connection  with 
the  Fokker-Planck  formalism.  The  model  encapsulates  the  coupled  dynamics  of  density  of  agent  population  with 
that  of  pheromone  laying,  following,  and  evaporation  in  a  2D  environment  and  carries  out  a  preliminary  numerical 
stability  analysis  of  the  solutions.  In  addition,  the  authors  also  present  a  distributed  control  law  in  continuous  time 
that  combines  gradient  following  for  pheromone  concentration  as  well  as  food  scent  with  random  motion  seen  in 
ants. 

The  second  scenario  relates  to  a  system  of  robots  interacting  via  non-linear  potentials  and  experiencing  random 
excitations.  The  proposed  framework  allows  calculation  of  probability  density  functions  of  states  of  the  robots, 
such  as  positions  and  velocities,  that  helps  analyze  the  stability  of  system  and  study  its  response  to  random 
excitations.  The  framework  allows  verification  of  the  notion  that  random  effects  in  swarm  systems  actually  help 
them  achieve  robust  self-organization  and  impart  adaptive  behaviors  in  changing  dynamic  conditions.  Although  this 
notion  has  existed  for  a  long  time,  a  comprehensive  mathematical  framework  to  study  the  effect  of  randomness  and 
to  actively  design  randomness  to  achieve  desirable  properties  has  been  absent  in  literature.  The  paper  formulates 
the  dynamics  of  two  robots  interacting  via  a  non-linear  potential  and  experiencing  random  forces.  The  notion  of 
adaptability  is  achieved  by  using  a  bistable  interaction  potential  where  the  two  stable  points  represent  the  two 
desirable  configurations  of  robots.  Adaptation,  in  this  case,  refers  to  the  ability  of  the  system  to  change  from  one 
configuration  to  the  other  due  to  change  in  environmental  conditions.  Bistable  interaction  potential  is  implemented 
by  introducing  a  quadratic  potential  (called  control  potential)  to  a  Morse  potential.  The  analytical  framework 
presented  in  this  paper  allows  the  study  of  the  effect  of  varying  noise  intensity  on  the  ability  of  the  system  to 
switch  from  one  stable  configuration  to  the  other  and  thus  presents  a  novel  approach  to  the  control  of  swarm 
systems  aimed  at  extracting  desirable  emergent  behavior  from  such  complex  systems. 


I.  Introduction 

Advances  in  communication,  computation,  and  sensing  systems  have  provided  unprecedented  capabili¬ 
ties  to  generate,  process,  and  communicate  vast  amounts  of  information.  From  the  perspective  of  military 
applications,  these  advances  have  the  potential  to  induce  a  paradigm  shift  in  the  military  Command 
and  Control  (C2)  by  providing  a  new  technological  framework  for  network  centric  warfare  capability. 
Highlights  of  the  framework  include  novel  mechanisms  for  sharing  information  and  enhanced  real-time 
situational  awareness  across  different  operational  groups,  resulting  in  increased  agility,  coordination, 
speed  of  command,  and  synchronization  in  C2.  From  the  standpoint  of  basic  research,  the  technological 
innovations  described  above  are  critically  dependent  on  fundamental  advances  in  the  area  of  complex 
networked  systems.  However,  as  affirmed  by  a  recent  National  Research  Council  report  by  Committee 
on  Network  Science  for  Future  Army  Applications  [2],  the  current  state  of  knowledge  in  this  field  lacks 
the  tools  to  that  can  accurately  predict  the  properties  of  large  scale  networked  systems.  In  particular, 
a  comprehensive  framework  comprising  rigorous  mathematical  and  computational  approaches  to  model, 
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predict,  and  control  the  behavior  of  such  systems  is  lacking  in  the  present  time.  In  this  context,  analytical 
approaches  using  ordinary  differential  equations  in  control-theoretic  frameworks  originally  developed  for 
simple  physical  systems  have  met  with  only  limited  success  in  the  case  of  complex  networked  systems 
primarily  because  of  two  fundamental  reasons:  (1)  the  inability  to  handle  the  higher  dimensionality 
that  characterizes  most  complex  systems,  and  (2)  the  incapacity  to  effectively  address  inherent  dynamic 
uncertainties.  Furthermore,  salient  characteristics  of  complex  networked  systems  include  their  ability  to 
adapt,  self-organize,  and  exhibit  other  emergent  behavior  which  often  fails  to  reveal  itself  in  an  analysis 
of  the  individual  sub-systems  using  traditional  methods.  In  summary,  complex  networked  systems  are 
certain  to  play  pervasive  roles  in  future  army  applications  and  the  lack  of  a  robust  analytical  framework 
to  model,  analyze,  and  control  their  unique  dynamical  behavior  represents  a  significant  challenge  from 
both  the  theoretical  and  applications  perspectives. 

Motivated  thus,  this  research  focuses  on  a  class  of  complex  networked  systems  namely  swarm  robotic 
systems  typically  represented  by  a  large  number  of  robots  in  mutual,  potentially  nonlinear,  interaction 
whilst  collectively  exhibiting  emergent,  group  behavior.  Swarm  systems  have  gained  considerable  attention 
among  the  robotics  community  primarily  due  to  three  reasons:  i)  robustness  of  these  systems  against 
failure  and  in  the  presence  of  uncertainty,  ii)  emergence  of  global  behavior  with  global  performance 
characteristics  that  cannot  be  achieved  by  the  resources  available  at  the  level  of  the  individual  unit,  and 
iii)  scalability  properties  due  to  the  unique  features  of  local  interactions  and  decision  making  processes. 
Major  challenges  in  developing  analytical  frameworks  include  nonlinearities  in  the  interactions,  typically 
high  dimensionality  of  the  state  space  of  the  full  system  comprised  of  several  subsystems,  and  uncertainties 
that  are  characteristic  of  the  sensing,  communication,  and  actuation  processes.  Additionally  it  is  well  known 
that  swarm  robotic  systems  exhibit  adaptive,  self-organizing  and  emergent  behavior.  However,  currently, 
the  ability  to  predict  such  global  behavior  from  the  analysis  of  local  subsystems  and  their  dynamical 
interaction  as  well  as  the  capability  to  control  the  emergent  behavior  to  achieve  desired  performance 
levels  are  conspicuously  absent. 

In  this  context,  paradigms  from  statistical  mechanics  which  deals  with  complex  systems  with  many 
degrees  of  freedom  and  several  states  offer  immense  potential.  Of  particular  note  is  the  Fokker-Planck 
formalism  that  aims  to  determine  probability  density  functions  for  the  system  evolution  and  hence  can 
be  used  to  determine  mean  values  of  dynamical  variables  rather  than  to  obtain  exact  trajectories.  In  the 
presence  of  uncertainty,  the  system  of  ordinary  differential  equations  representing  a  dynamical  system 
is  replaced  by  a  stochastic  differential  system,  the  simplest  example  of  which  is  the  Langevin  equation. 
Fundamental  results  in  stochastic  analysis  provide  a  canonical  method  of  obtaining  a  partial  differential 
equation  satisfied  by  the  probability  density  function  associated  with  the  Langevin  equation  (see,  for 
instance  [7]).  All  available  information  about  the  dynamical  response  of  a  system  driven  by  random 
excitation  is  encapsulated  in  the  probability  density  function  which  identically  satisfies  this  deterministic 
partial  differential  equation  -  the  Fokker-Planck  equation  ([7],  [26],  [29],  [18]).  Hence  solving  this  equation 
allows  complete  statistical  characterization  of  all  dynamical  random  variables  associated  with  the  system 
thereby  categorically  solving  the  problem.  While  the  analysis  of  dynamical  uncertainty  using  partial 
differential  equations  satisfied  by  probability  density  functions  is  a  classical  technique  that  has  found 
success  in  many  modem  applications  such  as  nanoelectromechanical  systems  [20],  [24],  the  Fokker- 
Planck  formalism  has  not  yet  been  comprehensively  investigated  in  the  context  of  swarm  robotics  systems 
and  this  serves  as  the  primary  motivation  for  the  approach  taken  in  this  research. 

In  broader  terms,  the  framework  adopted  in  this  report  has  immense  potential  to  open  up  various  exciting 
avenues  of  research  in  complex  networked  systems.  Of  particular  interest  are  phenomena  that  arise  uniquely 
due  to  the  synergetic  influence  of  noise  and  nonlinearity  and  are  being  actively  pursued  in  the  context  of 
dynamical  systems.  A  case  in  point  is  the  counterintuitive  phenomenon  of  stochastic  resonance  wherein 
the  addition  of  noise  to  certain  coupled  nonlinear  systems  leads  to  coherently  amplified  response  [6].  Yet 
another  instance  is  the  influence  of  stochastic  forcing  on  complex  systems  such  as  coupled  microresonator 
arrays  where  the  effects  of  noise  on  interesting  phenomena  such  as  energy  localizations  and  noise  induced 
enhanced  response  have  recently  been  studied  analytically  using  the  Fokker-Planck  approach  [23],  [21], 
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[22].  Whether  such  phenomena  arise  in  swarm  robotic  systems  and  if  they  do,  whether  they  could  be 
harnessed  to  practical  advantage  are  all  exciting  open  questions  that  can  be  pursued  within  the  proposed 
framework. 

The  rest  of  the  report  is  set  as  follows.  In  Section  2,  we  approach  the  dynamics  and  control  of  a 
robotic  swarm  motivated  by  an  analytical  model  of  ant  foraging  in  mathematical  biology.  A  continuous 
time  reaction-diffusion  model  obtained  from  a  system  of  Langevin  equations  that  essentially  yields  a 
Fokker-Planck  equation  for  the  robot  density  is  presented.  The  model,  which  yields  equations  similar 
to  the  Keller-Segel  model  of  chemotaxis,  encapsulates  the  coupled  dynamics  of  the  density  of  agent 
population  with  that  of  pheromone  laying,  following,  and  evaporation  in  two  spatial  dimensions.  Next, 
we  carry  out  a  preliminary  numerical  stability  analysis  of  the  solutions  obtained  from  the  model.  As  a 
highlight  of  our  treatment,  we  also  present  a  distributed  control  law  in  continuous  time  that  combines 
gradient  following  for  pheromone  concentration  as  well  as  food  scent  with  random  motion  seen  in  ants. 
In  addition,  we  apply  recent  analytical  results  on  stability  to  uniform  solutions  of  the  model.  In  Section 
3,  we  begin  by  applying  the  Fokker-Planck  formalism  to  a  system  of  robots  interacting  via  non-linear 
potentials  and  experiencing  random  excitations.  We  then  focus  on  the  dynamics  of  two  robots  interacting 
via  the  deterministic  Morse  potential  and  driven  by  white  noise.  The  notion  of  adaptability  is  achieved  by 
using  a  bistable  interaction  potential  where  the  two  stable  points  represent  the  two  desirable  configurations 
of  robots.  Adaptation,  in  this  case,  refers  to  the  ability  of  the  system  to  change  from  one  configuration 
to  the  other  due  to  change  in  environmental  conditions.  Bistable  interaction  potential  is  implemented  by 
introducing  a  quadratic  potential  (called  control  potential)  to  the  Morse  potential.  The  mathematical  tools 
developed  in  this  paper  allows  the  study  of  the  effect  of  varying  noise  intensity  on  the  ability  of  the 
system  to  switch  from  one  stable  configuration  to  the  other.  Concluding  remarks  are  collected  along  with 
a  discussion  of  directions  of  future  research  in  Section  4. 

II.  Control  Laws  for  a  Robotic  Swarm  Emulating  Ant  Foraging  Behavior 

A  significant  challenge  in  swarm  robotics  is  the  design  and  control  of  a  robotic  swarm  capable  of 
adaptive  behavior  dictated  by  local  communication  in  uncertain  environments.  Despite  the  problem  being 
a  formidable  one,  the  potential  pay  off  is  attractive  as  robotic  swarm  systems  play  a  fundamental  role  in 
a  broad  spectrum  of  cutting  edge  technological  applications  including  space  exploration,  nano-scale  drug 
delivery,  large-scale  surveillance,  and  search  and  rescue  operations.  In  all  these  applications  one  essentially 
deals  with  a  system  comprising  a  large  number  of  robots  working  towards  achieving  prescribed  objectives. 
The  dynamics  of  the  system  is  at  once  nonlinear  and  stochastic  due  to  the  inevitable  uncertainties  in  the 
environment  as  well  as  the  complexity  of  the  interactions  between  the  individual  elements.  The  success 
of  robotic  swarm  systems  as  a  technological  paradigm  is  ultimately  dependent  on  our  ability  to  extract 
desirable  adaptive  behavior  (such  as  self-organization)  from  such  complex  dynamical  systems  in  the 
absence  of  external  control  interference. 

Given  their  remarkable  propensity  to  adapt  to  rapidly  changing  environments,  the  dynamics  of  biological 
systems  provide  fundamental  insights  in  this  context.  Of  particular  interest  is  the  modeling  of  robotic  swarm 
inspired  by  the  behavior  of  ant  colonies  [13],  [11],  [17].  For  instance,  ants  are  capable  of  choosing  the 
best  nest  sites  out  of  multiple  possibilities  during  migration  and  reproductive  swarming.  Furthermore,  ant 
colonies  efficiently  establish  shortest  distance  trails  to  food  sources  that  minimize  foraging  time.  This  is 
achieved  by  the  deposition  of  a  chemical  known  as  pheromone  by  ants  returning  to  the  nest  from  profitable 
food  sources.  Successive  foragers  from  the  nest  are  motivated  to  follow  the  pheromone  trails  and  in  turn 
they  reinforce  the  trail  by  depositing  more  of  the  chemical  depending  on  the  continued  attractiveness  of 
the  food  source.  Various  analytical  models  have  been  proposed  to  study  the  social  behavior  of  biological 
systems  in  general  and  ant  colonies  in  particular  (see,  for  instance  [8],  [16]).  Furthermore,  the  constructive 
role  played  by  noise  in  collective  decision  making  during  foraging  has  been  reported  in  interesting  recent 
work  [4]. 

In  this  section,  we  discuss  the  analysis  and  control  of  a  swarm  of  mobile  robots  based  on  an  analytical 
model  [25]  inspired  by  ant  foraging  behavior.  Specifically,  we  present  a  distributed  control  law  in  continu- 
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ous  time  that  combines  gradient  following  for  pheromone  concentration  as  well  as  food  scent  with  random 
motion  seen  in  ants.  In  addition,  we  provide  a  continuous  time  model  for  pheromone  laying  in  a  two 
dimensional  environment.  Extensive  simulation  studies  confirm  emergent  behaviors  seen  in  ant  systems 
such  as  trail  formation  and  convergence  to  single  food  site.  Furthermore,  we  also  discuss  the  effect  of 
randomness  on  robustness  of  convergence  to  a  single  food  site.  It  is  noted  that,  in  contrast  to  previous  grid 
based  approaches,  a  continuous  robot  model  is  employed  here.  To  the  best  of  our  knowledge,  generating 
ant-like  behavior  in  a  swarm  of  mobile  robots  using  continuous  control  laws  has  not  been  reported  in  the 
literature  and  is  a  highlight  of  the  present  work. 


A.  Analytical  Model 

The  two  dimensional  dynamics  of  a  swarm  of  mobile  robots  dictated  by  pheromone  concentrations 
may  be  essentially  captured  in  a  two  component  reaction-diffusion  model  comprising  coupled  evolution 
equations  for  the  robot  density  a(r,t)  and  the  pheromone  density  b(r,t),  where  r  e  R2  represents 
position  in  two  dimensional  environment  and  t  represents  time.  The  model  described  here  is  a  modified 
version  of  the  one  first  discussed  in  [28].  The  robot  dynamics  evolves  on  a  two  dimensional  surface 
described  by  a  potential  U  (r,  t)  which  is  the  sum  of  contributions  from  the  environment  Uenv  and  the 
pheromone  deposit  Uph ■  Indeed  the  robots  are  able  to  alter  Uph  locally  by  the  production  of  pheromone 
and  the  explicit  functional  relationship  between  the  pheromone  density  and  the  potential  is  postulated  to  be 
Uph  =  —gb  (r,  t)  where  g  is  a  dimensional  constant.  The  environmental  potential  Uenv  encodes  information 
about  desirable  locations  in  the  space,  and  in  the  context  of  ant  foraging  can  be  considered  as  potential 
whose  gradient  will  drive  the  ants  to  the  food  sites.  Assuming  localized  pheromone  production  at  a  rate 
q  at  location  R, ,  allowing  for  the  pheromone  to  independently  diffuse  over  the  surface  with  a  diffusion 
coefficient  Db  as  well  as  to  evaporate  linearly  at  a  rate  7,  the  evolution  of  b  (r,  t)  may  be  modelled  as  a 
reaction-diffusion  equation  given  by: 

r)h  (r  /  i 

— — —  =  -7 b  (r,  t)  +  q^5(r  -  Ri  (t))  +  DbV2b  (r,  t) .  (1) 

i= 1 

where  R.t  represents  the  position  of  robot  i  (1  ranges  from  1  to  N ),  5  is  the  Dirac  delta  function  and  V2 
stands  for  the  Laplacian  in  cartesian  co-ordinates.  The  first  term  in  Eqn  (1)  represents  the  evaporation 
of  pheromone,  the  second  term  represents  deposition  of  pheromone  at  location  RL  and  the  third  term 
represents  diffusion  of  pheromone. 

The  stochastic  motion  of  robot  i  in  the  potential  U  (r,  t)  under  the  influence  of  white  noise  process 
dW  of  intensity  a  and  damping  coefficient  (3  may  be  represented  by  the  first  order  system  of  stochastic 
differential  equations  given  by: 

Ri  =  Vi) 

vi  =  -Pvi-VU(r,t)\Ri  +  adW. 

where  V  represents  the  gradient  with  respect  to  position  r.  In  the  Einstein-Smoluchowski  limit[28],  the 
inertia  term  in  Eqn. (2)  may  be  neglected  to  obtain: 

PRi  =  -VU(r,t)\Ri  +  adW.  (3) 


Substituting  the  explicit  expression  for  the  potential  U  (r,  t)  in  terms  of  b  (r,  t)  into  Eqn. (3),  and  considering 
Uenv  to  be  uniformly  distributed  (when  no  extra  information  regarding  the  environment  is  available),  one 
obtains: 


R,  = -jyjbir,t)\Rl+j,dW. 


(4) 
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Now,  the  discrete  description  of  the  individual  robots  is  replaced  by  a  continuum  description  characterized 
by  the  robot  density  a(r,t).  Accordingly,  the  pheromone  growth  term  (the  second  term  in  the  RHS  of 
Eqn(l))  is  replaced  by: 

N 

q'EHr-Mt))  qNp(r,t )  ->  qa(r,t),  (5) 

2=1 

where  p  (r,  t)  represents  the  probability  density  of  finding  a  robot  at  point  r  at  time  t.  Using  this  description, 
the  Fokker-Planck  equation  for  the  robot  density  a  (r,  t)  corresponding  to  the  stochastic  differential 
equation  (4)  can  be  derived  as: 

i)U  /jl' 1  =  ^ v  '  (-«  (A  t)  Vb  (r,  t))  +  DaX72a  (r,  t) ,  (6) 

where  Da  stands  for  the  diffusion  coefficient  corresponding  to  the  robot  density  a  (r,  t),  and  V-  represents 
the  divergence  operator.  Finally,  introducing  the  robot  density  (as  obtained  in  the  expression  (5))  into  the 
evolution  equation  for  the  pheromone  density  Eqn.  (1),  we  obtain: 

=  ~^b  (r>  t)  +  Qa  (r>  t)  +  DhSJ2b  (r,  t) .  (7) 

Equations  (6)  and  (7)  represent  the  coupled  evolution  equations  for  the  robot  density  and  the  pheromone 
density  and  as  such  form  the  analytical  basis  of  the  model.  It  may  be  noted  that  the  coupled  Equations  (6) 
and  (7)  represent  the  Keller-Segel  model  for  bacterial  chemotaxis  that  models  the  spatio-temporal  evolution 
of  population  of  bacteria  under  influence  of  chemicals  that  they  deposit  and  external  enviornmental  stimuli. 
The  model,  first  proposed  by  Keller  and  Segel  in  1970  [12],  has  been  researched  much  in  fields  of 
mathematics,  physics  and  mathematical  biology  [19],  [5],  [9],  [3],  [10],  [1],  [27],  [30].  However,  its  direct 
implementation  to  obtain  robot  control  laws,  to  the  best  of  authors’  knowledge,  has  been  absent. 

B.  Numerical  Stability  Analysis 

It  is  of  interest  to  examine  the  stability  of  solutions  to  the  coupled  system  of  equations  6  and  7  with 
respect  to  the  various  co-efficients.  Of  particular  significance  is  the  behavior  of  solutions  with  respect  to 
variations  in  magnitude  and  sign  of  g  as  well  as  changes  in  the  rate  of  pheromone  evaporation  7.  We 
note  that  positive  values  of  g  imply  an  attractive  pheromone  potential  and  will  result  in  a  highly  localized 
swarm.  On  the  other  hand,  a  negative  value  of  g  indicates  a  repulsive  pheromone  potential  profile.  The 
spatial  distribution  of  pheromone  obtained  from  numerical  simulations  of  Eqns.  6  and  7  at  the  end  of 
time  t  =  200  units  for  different  values  of  g  and  7  are  provided  in  Figures  1(a)- 1(d).  The  corresponding 
robot  distributions  are  provided  in  Figures  2(a)-2(d).  As  observed  from  Fig.  1(a),  for  the  case  of  negative 
values  of  g  and  no  evaporation  (7  =  0),  the  pheromone  density  shows  a  highly  ordered  pattern  of  spatial 
diffusion  at  the  end  of  the  simulation.  Correspondingly,  from  Fig. 2(a),  the  robots  show  a  high  degree  of 
spatial  dispersion.  However,  in  the  presence  of  a  relatively  small  evaporation  rate  (7  =  0.1)  the  robots 
begin  to  show  a  tendency  to  converge  (see  Fig. 2(b))  whereas  the  pheromone  is  much  more  even  distributed 
(Fig.  1(b)).  However,  for  positive  values  of  g,  which  may  be  understood  as  a  strongly  attracting  potential, 
and  relatively  low  evaporation  rate  (7  =  0.1),  the  pheromone  density  is  highly  concentrated  around  a 
point  (Fig.  1(c)).  Accordingly,  the  robots  cluster  around  this  point  (Fig. 2(c))  in  a  stable  configuration. 
Interestingly,  a  high  evaporation  rate  (7  =  10)  under  the  same  conditions  appears  to  induce  instability  in 
the  cluster  (see  Figs.  1(d)  and  2(d)).  This  can  be  explained  by  the  fact  that  due  to  high  evaporation  rate, 
the  effect  of  pheromone  on  dynamics  of  robots  is  lesser  pronounced  than  the  effect  of  randomness  (or 
diffusion)  on  the  dynamics.  Hence,  the  robots  appear  to  disperse.  These  numerical  results  indicate  that 
the  dynamical  stability  of  the  configurations  emerging  from  the  analytical  model  are  critically  dependent 
on  the  physical  parameters  of  the  system  of  governing  equations. 
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Fig.  1.  Spatial  pheromone  distribution  at  t  =  200  for:(a)  g=-5;7=0  (b)  g=-5;7=0.1  (c)  g=5;7=0.1  (d)  g=5;7=10 
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Fig.  2.  Spatial  robot  distribution  at  t  =  200  for:(a)  g=-5;7=0  (b)  g=-5;7=0.1  (c)  g=5;7=0.1  (d)  g=5;7=10 


C.  Stability  of  Uniform  Steady  State  Solutions 

As  mentioned  earlier,  stability  of  solutions  of  Keller-Segel  model  has  been  well  studied  in  literature 
especially  with  respect  to  blow  up  or  spike  formation  of  solutions.  Of  particular  interest  is  the  steady  state 
solutions  of  the  Keller-Segel  equations  that  result  into  a  uniform  density  of  chemicals  (pheromone)  and 
bacteria  (robots  or  ants).  There  has  been  some  recent  research  carried  out  on  this  aspect  that  have  aimed 
at  obtaining  realtionships  between  parameters  of  Keller  -  Segel  model  such  as  diffusion  coefficients  of 
pheromone  and  ants,  and  pheromone  evaporation  and  deposition  rate.  In  this  respect,  readers  are  referred 
to  papers  by  Payne  and  Straughn  [19]  and  Solis  et.  al.  [30]. 

A  preliminary  study  was  carried  out  to  verify  our  robot  control  model  using  the  results  published 
in  the  above  two  papers.  In  order  to  verify  the  results,  various  parameters  were  chosen  that  sastified 


Fig.  3.  Distribution  of  pheromone  (ID  case)  at:  (a)  t=0  (b)  t=67  (c)  t=133  (d)  t=200 
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Fig.  4.  Distribution  of  pheromone  (2D  case)  at:  (a)  t=0  (b)  t=67  (c)  t=133  (d)  t=200 


the  conditions  for  steady  state  uniform  solutions  reported  in  the  above  literature.  Figure  3  shows  the 
evolution  of  pheromone  density  in  one-dimensional  space  (y-axis  represents  pheromone  concentration 
and  x-axis  represents  the  ID  space)  for  disfferent  time  instants.  It  can  be  seen  that  as  time  elapsed  the 
spatial  distribution  of  pheromone  approached  towards  a  uniform  distribution.  The  results  are  as  well  valid 
for  the  2-dimensional  case  as  are  shown  by  Figures  4  and  5  that  respectively  represent  the  evolution 
of  pheromone  and  robots.  It  can  be  seen  that  robots  start  off  initially  from  a  clustered  position  and  get 
distributed  uniformly  as  time  elapses.  The  simulation  results  validate  the  derivation  of  individual  robot 
control  law  as  stochastic  differential  equation  given  by  (4)  corresponding  to  the  Fokker  Planck  equation 
given  by  (6). 

D.  Discrete  Ant  Behaviors 

While  the  two  component  reaction-diffusion  model  presented  earlier  provides  a  reasonable  analytical 
basis  for  ant  foraging,  in  order  to  better  reproduce  the  phenomenon  in  simulations,  we  need  to  consider 
some  of  the  discrete  behaviors  exhibited  by  ant  colonies  while  foraging.  Of  particular  interest  are  the 
following  two  types  of  behaviors:  i)  search  for  food;  and  ii)  transportation  of  food.  In  the  ’’search  for 
food”  behavior,  an  ant’s  motion  is  goverened  by  three  factors:  a)  environmental  potential  gradient  (VUenv) 


9 


Fig.  5.  Distribution  of  robots  (2D  case)  at:  (a)  t=0  (b)  t=67  (c)  t=133  (d)  t=200 

which  can  be  construed  as  a  gradient  generated  by  the  smell  of  food;  b)  gradient  of  pheromone  (VUph) 
laid  by  other  ants;  and  iii)  inherent  randomness  which  allows  the  ant  to  explore.  The  dynamics  of  motion 
of  ants  engaged  in  this  behavior  is  given  by  Equation  (4).  Furthermore,  the  intensity  of  randomness  is 
dependent  on  whether  an  ant  has  picked  any  gradient  on  either  environmental  potential  or  pheromone 
potential.  If  an  ant  has  not  picked  the  gradient,  then  it  will  have  certain  specific  amount  of  randomness  in 
its  motion.  This  randomness  allows  the  ant  to  explore  the  space  and  search  for  the  food  site  till  it  picks 
up  the  gradient.  This  kind  of  randomness  behavior  is  consistent  with  behavior  seen  in  actual  ants. 

During  the  ’’transportation  of  food”  behavior,  the  ant  carries  the  food  back  to  the  nest  while  laying  the 
pheromone  on  the  path.  In  accordance  with  actual  ant  system  behaviors,  the  ants  are  assumed  to  move 
with  a  constant  velocity  in  straight  line  from  the  food  site  to  the  nest. 

E.  Simulation  Results:  The  Base  Case 

In  this  section  we  present  the  simulation  results  for  the  base  case  of  very  low  noise  intensity.  Randomness 
is  introduced  into  the  dynamics  by  passing  the  effective  potential  U  (r,  t)  through  a  white  Gaussian  noise 
channel  using  the  MATLAB  function  ’awgn’.  This  is  equivalent  to  adding  white  noise  to  the  effective 
potential.  The  noise  intensity  is  modulated  by  specifying  the  signal-to-noise  ratio  (SNR)  of  the  channel. 
For  the  base  case,  a  relatively  high  value  of  SNR  =  10776  is  considered  and  therefore  the  dynamics  is 
almost  deterministic  in  this  case.  The  simulation  is  run  for  time  t  =  200  units  with  a  time  step  of  0.04. 
The  nest  is  represented  by  the  red  star  (with  grid  co-ordinates  (55,20))  whilst  the  blue  stars  indicate 
the  two  food  sources  with  grid  co-ordinates  (30,50)  and  (70,50)  respectively.  Uenv  is  considered  to  be 
unifrom  for  the  simulations.  The  food  source  located  at  (70,  50)  is  considered  to  be  superior  because  of 
its  proximity  to  the  nest.  Figures  6(a)-6(d)  represent  snapshots  of  the  traffic  and  time  t  =  0,  t  =  Q)  200, 
t  =  (|)  200  and  t  =  200  respectively.  The  formation  of  trails  is  quite  clear  from  the  simulation  results. 
In  addition,  snapshots  of  the  pheromone  density  in  both  the  trails  at  each  of  the  time  instances  indicated 
above  are  also  furnished  in  Figures  7(a)-  7(d)  using  the  color  code  indicated  in  the  figures.  It  is  interesting 
to  note  the  disparity  in  traffic  intensity  between  the  food  sources  with  the  superior  food  source  being  the 
clear  favourite.  This  is  further  demonstrated  in  Figure  8  wherein  the  number  of  robots  visiting  each  of 
the  food  sources  is  plotted  as  a  function  of  time. 

F.  Simulation  Results:  SNR=2 

As  a  step  towards  studying  the  effects  of  noise,  in  this  section  we  present  simulation  results  for  the  case 
of  high  intensity  noise  (SNR=2)  influencing  the  dynamics  once  a  robot  has  picked  up  a  gradient.  Similar 
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Fig.  7.  Trails  to  the  food  sources  as  demonstrated  by  pheromone  density  (Base  Case:  SNR  =  10E6)  at:(a)  t=0  (b)  t=67  (c)  t=1 33 
(d)  t=200 


to  the  base  case,  Figures  9(a)  to  9(d)  represent  snapshots  of  the  traffic  and  time  t  =  0,  t  =  (|)  200, 
t  =  (|)  200  and  t  =  200  respectively.  The  corresponding  trails  are  depicted  in  Figures  10(a)- 10(d)  and 
the  number  of  robots  visiting  the  two  food  sources  as  functions  of  time  are  plotted  in  Figure  1 1 . 

G.  Effects  of  Varying  Noise  Intensity 

To  better  understand  the  effect  of  varying  the  intensity  of  noise  on  the  traffic,  in  this  section  we  present 
results  for  the  time  variation  of  the  number  of  robots  visiting  the  two  food  sources  for  the  cases  of  SNR=3, 
SNR=10,  SNR=100  and  SNR=1000  in  Figures  12(a),  12(b),  12(c)  and  12(d)  respectively.  Comparing  the 
base  case  (Figure  8)  with  the  cases  of  higher  noise  intensity  (Figures  11  and  12)  we  first  note  that 
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Fig.  8.  Number  of  robots  visiting  Food  Source  1  (red)  and  Food  Source  2  (blue). 


(c)  (d) 

Fig.  9.  Configuration  of  Robots  for  SNR=2  at:(a)  t=0  (b)  t=67  (c)  t=1 33  (d)  t=200 


addition  of  noise  does  not  fundamentally  alter  the  preference  of  the  better  food  source  with  food  source  1 
continuing  to  be  favored.  However,  in  the  base  case  of  very  low  intensity  noise  the  traffic  to  food  source 
1  shows  an  abrupt  increase  at  around  t  =  25  while  an  increase  in  traffic  towards  food  source  2,  albeit  to 
a  lesser  extent,  occurs  between  t  =  40  and  t  =  60.  One  observes  a  qualitative  difference  in  this  feature  as 
the  noise  intensity  is  increased.  As  can  be  seen  from  Figures  11  and  12,  at  higher  intensity  of  noise,  the 
traffic  to  both  the  sources  picks  up  at  much  earlier  times  and  also  shows  growth  almost  proportional  to 
the  time  elapsed. 


III.  Fokker-Planck  Equation  for  Morse  Potential 

In  this  section,  we  present  the  Fokker-Planck  formalism  within  which  a  swarm  robotic  system  can 
be  modeled  and  analyzed.  The  swarm,  comprising  a  system  of  N  coupled  robots  undergoing  collective 
motion  in  which  robot  i  experiences  a  total  interaction  potential  U  (V/,;)  involving  attractive  and  repulsive 
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(C)  (d) 

Fig.  10.  Trails  to  the  food  sources  as  demonstrated  by  the  pheromone  density  for  SNR=2  at:(a)  t=0  (b)  t=67  (c)  t=1 33  (d)  t=200 


Fig.  11.  Number  of  robots  visiting  Food  Source  1  (red)  and  Food  Source  2  (blue)  for  case  SNR=2. 


forces  resulting  from  pairwise  potential  U  (qtJ ) ,  a  control  input  dissipative  force  f  (p^  and  white  noise 
process  dW  of  strength  a,  may  be  represented  as: 

ft  =  Pu  Pi  =  Ui  =  f  ( Pi )  -  ViU  (ft)  +  crdWi .  (8) 


Artificial  potential  functions  have  been  used  frequently  [15],  [14]  for  modeling  of  intractions  in  swarm 
robotic  systems.  Compare  the  above  system  of  equations  (8)  with  a  standard  Ito  SDE,  dx  =  Adt  +  BdW, 
where  x  is  the  state  vector,  A  is  the  drift  vector  and  B  is  the  diffusion  matrix.  The  Fokker-Planck  equation 
for  the  probability  density  Pj  governing  the  dynamics  of  robot  i  may  be  derived  as: 


m 

dt 


d 


dx , 


[AjPi]  +  ^ 


d 2 


j,kdxjdxk 


(9) 


Equation  (9)  provides  a  probabilistic  description  of  the  dynamics  of  each  element  of  the  swarm.  One  of  the 
objectives  of  this  research  to  study  the  influence  of  noise  on  emergent  behavior  of  swarm  robotic  systems. 
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Fig.  12.  Number  of  robots  visiting  Food  Source  1  (red)  and  Food  Source  2  (blue)  for:  (a)  SNFt=3  (b)  SNFt=10  (c)  SNFt=100  (d) 
SNFM000 


Contrary  to  common  intuition,  noise  can  sometimes  play  a  constructive  role  in  nonlinear  systems.  In 
particular,  as  evidenced  by  the  phenomenon  of  stochastic  resonance,  the  effects  of  noise  may  be  harnessed 
to  advantage  in  applications  [6].  In  interesting  recent  work,  Meyer  and  co-workers  have  suggested  that 
hitherto  unobserved  adaptive  behavior  in  ant  colonies  can  be  explained  using  noise  induced  transitions 
between  stable  equilibrium  points  of  a  double-well  potential  [4]. 

Inspired  by  the  above  biological  example,  we  carried  out  preliminary  research  to  study  the  noise  induced 
transitions  in  the  robotic  swarm  using  a  two-robot  system.  The  central  idea  is  to  introduce  a  control 
potential  to  the  regular  pairwise  interaction  potential  resulting  in  an  altered  effective  potential.  Altering  the 
resultant  interaction  potentials  using  control  interventions,  as  proposed,  results  in  new  equilibrium  points 
for  the  resultant  effective  potential.  Concomitantly,  the  stationary  densities  obtained  from  the  Fokker- 
Planck  equations  will  have  peaks  around  these  equilibrium  points.  The  control  objective  of  extracting 
desirable  adaptive  behavior  will  be  realized  by  driving  the  system  to  these  equilibrium  points  in  the 
probabilistic  sense. 

As  a  proof  of  concept  for  the  ideas  introduced  above,  we  now  consider  a  two  robot  system,  wherein 
the  mutual  interaction  is  represented  by  the  Morse  potential  given  by  the  following  equation: 

Uint  =  De(  l-e-^-^)2  (10) 

A  harmonic  oscillator  potential  as  the  control  intervention  is  given  by: 

Ucon  =  be{r  -  rc)2  (11) 

where  r  represents  the  inter-robot  distance.  The  resultant  effective  potential  is  given  by  U  =  Uint  +  Ucon. 
The  Morse  potential  is  a  well  known  potential  for  modelling  inter-robot  interaction  and  Figure  13  shows  a 
plot  of  Morse  potential  with  a  minimum  at  inter-robot  distance  of  three  (3)  units.  Addition  of  the  harmonic 
oscillator  potential  results  in  an  extra  ’’well”  as  can  be  seen  in  Figure  14  at  inter-robot  distance  of  9.5. 
Hence,  the  resultant  potential  has  two  minima:  one  at  inter-robot  distance  of  3  and  another  at  inter-robot 
distance  of  9.5. 
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Interaction  Potential 


Fig.  13.  Morse  potential  plotted  against  inter-robot  distance 


Interaction  Potential 


Fig.  14.  Resultant  effective  potential  plotted  against  inter-robot  distance 


Numerical  simulation  of  the  two-robot  system  undergoing  the  dynamics  given  by  Equation  8  under  the 
influence  of  combined  potential  and  noise,  as  demonstrated  in  Figure  15,  clearly  shows  the  transition  of  the 
system  between  equilibrium  states  for  a  noise  intensity  of  25.  The  key  points  to  be  noted  here  are:  (1)  new 
equilibria,  corresponding  to  desirable  emergent  behavior,  can  be  introduced  using  control  potentials  and 
(2)  optimal  levels  of  noise  intensity  may  be  used  to  control  transitions  of  the  system  between  equilibrium 
states. 

Finally,  we  derive  the  Fokker  -  Planck  equation  corresponding  to  the  stochastic  differential  equation 
governing  the  dynamics  of  the  two-robot  system  under  the  influence  of  an  effective  potential  given  by  the 
sum  of  Morse  potential  and  harmonic  oscillator  potential  as: 

=  -y%  +  fa + 2“B.  ('■“(r_")  -  '■2“(r_"))  +  *-]  %  + cP  +  y  <12> 

Here,  y  —  r,  and  P  represents  the  probability  density  of  r  and  y,  P(r.  y).  Steady  state  solutions  to 
the  above  Fokker- Planck  equation  may  be  utilized  to  analytically  characterize  the  transitions  between 
equilbria  that  were  numerically  demonstrated  earlier  in  this  section.  On  the  other  hand,  while  Fokker- 
Planck  equations  for  nonlinear  systems  are  notoriously  hard  to  solve,  approximate  techniques  such  as 
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Time 


Fig.  15.  Inter-agent  distance  plotted  against  time 


the  method  of  moment  approximations  or  cumulant  expansions  may  be  used  to  obtain  the  time  evolution 
of  averaged  dynamical  quantities  from  the  Fokker-Planck  equation.  These  investigations  are  currently 
underway  and  will  be  reported  in  future  publications. 

IV.  Conclusions  and  Future  Work 

In  Section  2  of  this  report,  we  presented  continuous  time  control  laws  for  robots  based  upon  the  Keller- 
Segel  model  of  bacterial  chemotaxis  designed  to  emulate  ant  foraging  behavior.  The  motion  control  law 
was  based  on  gradient  following  that  attempted  to  minimize  an  effective  potential  which  is  the  sum  of 
an  environmental  potential  and  the  potential  arising  due  to  pheromones.  Furthermore,  the  control  law 
incorporated  a  random  component  that  allowed  the  robots  to  stochastically  explore  the  spatial  domain  for 
food  sources.  A  continuous  time  model  of  pheromone  dynamics  that  includes  pheromone  deposition, 
diffusion,  and  evaporation  was  also  presented.  The  discrete  control  law  formulated  on  the  basis  of 
correspondence  of  Fokker-Planck  equation  to  stochastic  differential  equation  was  validated  based  upon 
solutions  of  Keller-Segel  model  reported  in  literature.  Finally,  a  set  of  discrete  behaviors  inspired  by 
foraging  ant  colonies  was  incorporated.  Exhaustive  simulations  demonstrate  two  of  the  emergent  behaviors 
seen  in  ant  foraging  :  i)  formation  of  trails,  and  ii)  increased  traffic  along  the  trail  that  leads  to  the 
superior  food  source.  In  addition,  preliminary  results  on  the  influence  of  noise  intensity  on  the  dynamics 
were  also  presented.  The  simulation  results  suggest  that,  interestingly,  higher  intensity  noise  encourages 
traffic  towards  the  better  food  source  in  a  manner  almost  directly  proportional  to  the  time  elapsed.  Several 
interesting  directions  of  future  research  emerge  from  these  results.  For  instance,  generalizations  of  the 
Keller-Segel  model  that  take  into  account  the  qualitative  differences  between  the  dynamics  exhibited  by 
foraging  ants  and  food  carrying  ants  are  important  topics  for  further  study.  On  the  other  hand,  while 
there  exists  a  vast  amount  of  mathematical  literature  pertaining  to  singular  behavior  of  solutions  to  the 
Keller-Segel  model  and  stability  of  solutions,  a  detailed  study  of  the  implications  of  the  results  for  the 
design  and  control  of  swarm  robotic  systems  is  absent  at  the  current  time.  The  initial  steps  in  this  direction 
taken  by  the  present  report  need  to  be  pursued  in  further  detail.  Important  open  questions  on  the  influence 
of  noise  intensity  on  swarm  dynamics  also  need  to  be  addressed  comprehensively. 

In  Section  3  of  this  article,  we  presented  a  direct  approach  to  studying  the  dynamics  of  a  collection  of 
interacting  entities  starting  with  the  Langevin  equations  and  then  deriving  the  Fokker-Planck  equation.  The 
novel  idea  of  extracting  desirable  emergent  behavior  by  altering  the  effective  potential  via  the  introduction 
of  control  potentials  was  introduced  and  applied  to  the  case  of  a  two-agent  system.  It  was  shown  that 
a  control  potential  of  the  harmonic  oscillator  type,  when  added  to  a  Morse  potential,  resulted  in  a  new 
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equilibrium  point  for  the  system.  Furthermore,  transition  between  the  two  equilibrium  points  for  an  optimal 
level  of  noise  intensity  was  shown  numerically.  The  Fokker-Planck  equation  for  this  scenario  was  derived 
and  an  analytical  characterization  of  the  state  transitons  is  currently  under  investigation.  In  addition,  the 
influence  of  noise  intensity  on  the  dynamics  is  being  studied  using  the  method  of  moment  approximations, 
the  starting  point  of  which  is  the  Fokker-Planck  equation  derived  above.  A  central  issue  in  the  use  of 
discrete  particle  models  in  swarm  robotics  is  the  smooth  passage  from  such  models  to  continuum  models. 
It  is  of  great  interest  to  investigate  the  effect  of  control  potentials  on  emergent  behavior  in  the  context  of 
continuum  models.  These  and  allied  questions  are  the  subject  of  ongoing  research  and  will  be  reported 
in  future  publications. 
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